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1.  Introduction 

In  this  paper  we  considers  the  numerical  solution  of  elliptic  par- 
tial differential  equations  in  spherical  domains.  There  are  numerous 
applications  In  engineering  and  the  physical  sciences  In  which  the  sol- 
ution of  a spherical  elliptic  equation  is  desired;  see  the  references 
In  (12}.  'When  all  the  functions  involved  are  spherically  symmetric 
(that  is,  they  depend  only  on  distance  from  the  center  of  the  domain), 
the  problem  can  be  replaced  by  an  equivalent  two-point  boundary  value 
problem.  The  resulting  problem  Is  singular,  but  nevertheless  has  a 
smooth  solution.  It  should  therefore  be  possible  to  approximate  the 
solution  accurately  usitig  the  Raylelgh-Rl tz  Calerkln  method  with  a 
piecewise  polynomial  subspace  on  a quaslunlform  mesh.  We  will  obtain, 
optimal-order  error  bounds',  showing  that  this  procedure  Is  theoretlci^ly 
well-founded.  Instead  of  the  usual  Sobolev  norms,  we^  use  norms.,  which 
are  appropriate  to  the  original  n-dimenslonal  setting  of  the  problem. 

The  problem  has  been  considered.  In  2 and  3 dimensions,  by  Russell 
and  Shamplne  (12].  They  obtain  error  bounds  for  approximation  proce- 
dures specially  designed  to  deal  with  the  apparent  singularity  at  the 
origin.  In  particular,  they  treat  collocation  In  which  the  basis  Is 
augmented  by  singular  basis  functions,  singular  patch  bases  (L-spllnes), 
and  a finite  difference  scheme  of  Jamet  (8)  designed  to  handle  the  sin- 
gularity. Crouzelx  and  Thomas  (11  and  Reddien  [11]  consider  this  prob- 
lem as  part  of  a wider  class  and  obtain  similar  results. 

Dupont  and  Wahlbin  [4]  and  Jesperson  (9)  have  analyzed  an  approxi- 
mation procedure  similar  to  ours,  and  have  obtained  error  bounds  of  op- 
timal order  in  the  usual  Sobolev  norms.  The  import  of  their  results, 
together  with  those  of  this  paper.  Is  that  no  special  measures  are  re- 
quired for  this  problem:  the  Raylelgh-Ritz  Calerkln  method  using 
high-order  piecewise  polynomial  spaces  on  a uniform  mesh  is  a highly 
effective  numerical  method. 

I 

[.  2.  Spherically  Symmetric  Elliptic  Problems 
! Let  B(a)  be  the  open  sphere  of  radius  a in  r",  l.e.. 


B(a)  = (x  e r"  | r(x)  < a}. 
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defined  on  B is  spherically  symmetric  If  It  depends  only  on  distance  r 
from  the  origin.  If  F and  Q are  spherically  symmetric  then  an  obvious 
symmetry  argument  shows  that  U Is,  too  (a  change  In  coordinate  systems 
by  rotation  around  any  axis  passing  through  the  origin  leaves  the  prob- 
lem, and  hence  Its  solution,  unchanged). 

Let  u(r),  f(r),  and  q(r)  be  functions  such  that 

F(jc)  - f(r(x)), 

Q(x)  - q(r(x)), 

U(x)  - u(r(x)). 

Then  u(r)  can  be  obtained  as  the  solution  of  the  singular  two-point 
boundary  value  problem 

- D(r"  ^Du)  + r''”^q(r)u  - r"  ^f(r),  0 < r < 1,  (3) 

Du(0)  - u(l)  - 0,  (4) 

%rttere  D = 

dr 

We  note  that  In  case  n ■ 3,  the  well-kno%m  change  of  variables 
V - ru  results  In  the  nonsingular  problem 

2 

-Dv+qv  - rf,  0<r<l, 
v(0)  - v(l)  - 0. 

However,  we  know  of  no  such  trick  In  two  dimensions! 

For  real-valued  functions  f,g  defined  on  (0,a)  we  let 


and 


a 

= r"  f(r)g(r)  dr, 
0 


Ilf  II 


B(a) 


(f.f) 


1/2 

B(a)* 


Furthermore,  we  define  j"(a)  (respectively  3^(3))  to  be  the  closure  of 

the  C functions  (respectively  the  C*  functions  which  vanish  In  a 
neighborhood  of  a)  with  respect  to  the  norm 


Ilf  II 


m,B(a) 


(a 


1/2 


u 


jk 
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A 
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Ue  assume  that  the  coefficient  q c C(I)  Is  such  that  there  exist 
positive  constants  1 and  A such  that 

for  all  u e J^(a),  where 

a 

lu,v]  . . = / r iDuDv  + quv)  dr. 

Bta;  Q 

Our  approximation-theoretic  results  In  the  spaces  J will  rely  on 
the  following  basic  fact,  proved  In  Courant  and  Hilbert  [2]. 

Lemma  1:  If  f c Jq(h) . then 

“'llB(a)  i • “““bo-  <‘> 

Inequalities  (5)  and  (6)  show  that  the  bilinear  form  [ , ]_  Is 

1 0 ® 

positive  definite  over  the  space  J>.  Thus,  for  each  f c J there  exists 
1 ” 

a unique  u c Jq,  called  the  generalized  solution  of  (3)  - (4),  such  that 
[u,v]g  » (f,v)g,  for  all  v e 

This  differential  equation  Is  exceedingly  well-behaved.  In  fact, 
we  know  that  the  size  of  the  solution  Is  bounded  In  terms  of  the  size  of 
the  data. 

Lemma  2:  There  exists  a constant  F such  that,  for  all  f e J®,  the  gen- 
eralized solution  u of  (3)  - (4)  satisfies 

II  D^u  II  g < r II  f II  g. 

Proof ; See  [13]. 

We  now  state  a variant  of  the  Sobolev  lemma.  Let  | x I denote  the 
greatest  integer  not  exceeding  x.  “ 

Lemma  3:  Let  u e J*”(a).  There  exists  a positive  constant  C such  chat. 


for  all  r e [0,a] , 
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Proof ; See  Friedman  [7J. 

In  the  cases  n - 2 or  3,  (7)  applies  with  m - 2. 

Let  S be  a f Inlte-dlmenslonal  subspace  of  ji.  The  function  a e S 
" On 

Is  the  Raylelgh-Rltz  Galerkln  (RRG)  approximation  to  u If 

)o»  ‘ll  v c S . 

no  n B n n 

Since  the  RRG  approximation  Is  the  projection  of  u on  S^  %flth  respect  to 
the  lnp'’r  product  we  have  the  error  bounds 


(u  - u,u  - u], 


Inf  lu  - V ,u  - V 1 , 

V cS  " " ® 

n n 


and  by  ( 5 ) , 


D(u-u)||_  < X"  A Inf  ||D(u-v)|l  . 

® V cS  " ® 


Let  A be  a partition  of  [0,1]: 

a:  0 - Xq  < x^  < X2  < ...  < - I. 

and  h = max  (x.  - x ).  We  assume  that  the  partition  Is  quaslunl- 
K1<N  ^ 

form,  meaning  that  the  global  mesh  ratio. 


*1  ■ *1-1 
M = max  — 

l<l,j<N  *j  “ *J-1 
Is  bounded  Independent  of  N. 

As  approximating  subspaces  we  will  use  spaces  of  plece%rlse  poly- 
nomials with  respect  to  the  partition  A: 

Sq(A,v)  = {s  e C'’[0,1])  I s(l)  - 0,  and 

s Is  a polynomial  of  degree  < k on  each  of 


the  Intervals  (x. 


_1*  1 i i £ N >. 
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3.  Spherical  Sp 1 Ine  Approx laatlon 

We  first  consider  approximation  by  polynomials  In  a neighborhood  of 
the  origin. 

n 

Lemma  4:  Let  v c J (a),  m > 1.  There  exists  a polynomial  T^v  of  order 
m satisfying 

for  all  0 _<  J < m. 


Proof ; Let  T^v  be  the  first  m terms  of  the  Taylor  series  for  v at  a, 
1 . e • . 

T v(x)  - v(a)  + Dv(a)(x-a)  + . . . + (x  - a)*"^ 

a m— 1 ! 

Since  - T v)  c J^Ca),  we  may  apply  Inequality  (6),  obtaining 

a u 


IId”‘^v  - T^v) 


B(a)  - 


V>  "B(a) 


- « "°"^"B(a)- 

To  prove  (9),  we  note  that  D^(v  - T v) (a)  - 0 for  all  0 < J < i»-l.  We 

a ~ 

may  therefore  apply  (6)  to  each  of  the  remaining  derivatives^  and  use 
the  result  for  to  obtain  the  result  for  D'^. 


Q.  E.  D. 

Theorem  For  each  integer  2 £ m £ k,  there  exists  a positive  constant 
K > K(k,m,n) , depending  on  the  global  mesh  ratio,  such  that  for  each 

f € J™  n Jq  there  exists  an  approximation  f c Sq(A,v)  to  f satisfying 

l|D^(f-f)||g  < K h"”J  IlD^fllg.  (10) 


Proof:  The  details  of  the  proof  may  be  found  In  [13].  We  will  sketch 
the  main  ideas  here. 


The  approximation  whose  existence  is  asserted  can  be  explicitly 

constructed.  Let  be  the  B-spllne  basis  functions  for  the  space 

k 1 

S^(a,v).  In  [3],  de  Boor  and  Fix  construct  a set  of  linear  functionals 
{ Ij)  dual  to  the  B-spllne  basis,  l.e.. 


X,(Bj) 


i - j 
1 y J * 
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The  approximation 


Ff  s £ X rf)B, 

^ 1-1  ‘ ‘ 

Is  called  Che  quasllnterpolant  of  f. 

The  functional  Is  a linear  combination  of  derivatives,  eval- 

uated at  some  point  c supporC(B^).  Moreover,  difference  approxima- 
tions may  be  used  Instead  of  derivatives. 

We  assume  Chat  for  all  1 such  chat  x^  Is  In  Che  support  of 

and  that  the  points  used  In  Che  corresponding  difference  approxima- 
tions are  all  contained  In  the  Interval  (0,Xj).  It  can  be  shown  that 
there  exists  a constant  C such  that,  for  all  such  1, 


X^(f)|  < C Ilf 


L (0,Xj) 


To  obtain  error  bounds  for  Che  quasllnterpolant  In  Che  weighted 

norm  ||  ||  , we  use  the  fact  Chat  Che  quasllnterpolant  of  any  polynomial 

B 

of  degree  < k Is  Chat  polynomial  [3].  Considering  first  Che  Interval 

(0,  x^. 


'a''  “b(,j)  i 


Let  R-f-Tf.  By  (11), 
1 


D^F  (T  f - f) 
A Xj 


B(Xj)‘ 


d-'f  (T  f - f) 


B(Xj)  - 


< £ |X,(R)|  II D-*  B,  II 


1 "B(Xj) 


< C ||R 


L“(0.xi)  «l«B(x^)* 


All  but  k of  the  basis  functions  vanish  In  (0,Xj).  For  those  that 
don't,  It  Is  readily  shown  that 

Together  with  the  Sobolev  lemma  and  the  bounds  for  R given  In  Lemma  4, 
this  yields 
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IID^F  (T  f - f) 
a a 


B(Xj)  ' 


K x: 


l|D"f 


B(Xj)' 


(12) 


A similar  argument  gives  the  same  result  for  the  Intervals 
^*l’*l+l^’  ^ ^ remaining  Intervals,  the  support  of  any 

basis  function  which  Is  nonzero  In  the  Interval  Is  bounded  away  from 
zero.  Using  this  fact,  an  approximation  result  analogous  to  (12)  la 
easily  obtained.  Combining  these  results  yields  (10). 


Q.  E.  D. 


i.  Error  Bounds  for  the  RRG  Approx Imatlon 

We  now  consider  the  error  In  the  RRG  approximation  to  the  general- 
ized solution  u of  (3)  - (4),  and  show  that  the  RRG  approximation 

ii  c Sj.(i,v)  approximates  u to  optimal  order. 

Theorem  2:  Let  u be  the  generalized  solution  of  (3) 
u e Sq'A.v)  be  the  RRG  approximation  to  u.  If  u c J™  J^, 
then 


l|D(u  - u)  11  g 

< 

1 

V| 

hm-l 

II  d"u  II  g 

(13) 

11  u - u II  g < 

(AK)^r  h™ 

II  0% 

"b- 

(14) 

Proof : The  error  bound  (13)  follows  immediately  from  (8)  and  the  result 
of  Theorem  1.  Inequality  (14)  also  follows,  via  Nltsche's  Crick,  from 
Theorem  1 [10]. 


- (4).  Let 
2 £ m < k. 


5.  Computational  Aspects  of  the  Method 

We  have  shown  that  the  Raylelgh-Ritz  approximation  u to  the  solu- 
tion u of  (3)  - (4)  Is  optimally  accurate  In  the  natural  norms  for  this 
problem.  We  would  like  to  think  of  this  approximation  as  inducing  a 
smooth,  accurate  approximation  of  the  solution  U of  the  original  problem 
(1)  - (2).  For  this  to  happen,  the  odd  derivatives  of  u must  vanish  at 
0.  Unless  we  Impose  this  requirement  on  the  subspace,  it  will  not  be 
f u 1 r i 1 1 ed  . 

There  are  two  ways  this  can  be  done.  We  can  simply  force  the  ele- 
nent«  of  the  space  Sp(A,v)  to  have  odd  derivatives  which  vanish  at  0. 
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Alternatively  we  consider,  rather  than  (3)  - (4),  the  two-point 
boundary  value  problem 

- D(|x|''"^Du)  + |x|"”^q(x)u  - |x|"~^f(x), 

-1  < X < I,  u(-l)  - u(l)  - 0, 

also  derived  from  (1)  - (2).  We  then  define  SS.(A,v)  to  be  the  space  of 

V 

C piecewise  polynomials  of  degree  < k (vanishing  at  -1  and  1)  with 
respect  to  a partition  which  Is  symmetric  about  0: 


- 0 < X,  < ...  « x^  • 1. 


where  x_^  " x^,  1 £ f £ N. 


It  can  be  shown  that  the  results  of  Theorems  1 and  2 hold  for  the  space 
SSq(4,v)  (131. 


The  B-spllne  basis  functions  then  have  a symmetry  property,  derived 
from  the  symmetry  of  the  mesh:  If  Bq,  B^  are  the  basis  functions 

(d+1  « dlm(SSQ(A,v) ) numbered  In  the  natural  lef t-to-rlght  order,  then 

B.  (-x)  - B.  .(x),  -1  < X < 1 , 0 < 1 < d. 

1 d— 1 “ 

d 

Clearly,  5 • E a B Is  even  If  the  vector  a Is  symmetric  about  Its 
1-0  ^ ‘ 

middle: 


a 


1 


a.  ^ , 0 < 1 < d. 

d— i — 


(15) 


The  symmetry  property  (15)  will  hold  for  the  RRG  approximation's 
coefficients,  even  though  we  do  not  Impose  It.  The  coefficients  are 
obtained  as  the  solution  of  the  linear  system 

AO  - f^,  (16) 

where  A - ♦ X “ 

1 . 

a^j  i / |x|""VdB^DBj  + qB^B^)  dx, 

1 , 

|x|””  t B.  dx. 

^ -1  ^ 

Clearly,  because  of  the  symmetry  of  the  data  and  the  basis  functions, 
the  matrix  A will  be  symmetric  about  the  alternate  diagonal  and  the 


ail 
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vector  vlll  be  symmetric  about  Its  middle.  This  shows  that  the  coef- 
ficients of  u will  satisfy  (15).  Thus  Q Is  even,  and  all  Its  odd  deri- 
vatives of  order  £ v vanish  at  0. 

It  does  not  cost  any  more.  In  work  and  storage,  to  use  the  (-1,1) 
pr-iblen  than  the  (0,1)  problem.  Because  of  the  symmetries  of  A,  only 
1/4  of  its  elements  need  to  be  computed.  Moreover,  using  an  algorithm 
of  Evans  and  Hatzopoulos  [5]  which  takes  advantage  of  symmetry  about  the 
alternate  diagonal,  the  equations  (16)  can  be  solved  in  half  the  time 
required  by  the  usual  band  Cholesky  algorithm. 

The  effect  of  numerical  quadrature  (used  to  compute  the  matrix  A 
and  the  right-hand  side  vector  £)  on  the  accuracy  of  the  RRC  approxima- 
tion has  been  analysed  by  Fix  for  nonsingular  problems  [6].  He  showed 
that  If  the  Integrals  are  computed  using  composite  Gaussian  quadrature 
with  k-1  points  In  each  Interval,  then  the  error  due  to  the  quadrature 
Is  asymptotically  as  small  as  the  discretization  error.  We  have  been 
able  to  show  that  k points  Is  sufficient  for  this  type  of  singular 
problem  [13],  but  conjecture  that  this  result  can  be  Improved,  and  that 
k-l  points  also  suffice  here.  The  numerical  results  of  the  next  section 
strongly  support  this  viewpoint. 


h.  Numerical  Results 

In  this  section  we  present  the  results  of  a numerical  experiment, 
which  Illustrates  the  utility  of  the  computational  procedure  analyzed  In 
the  previous  sections.  Following  Russell  and  Shampine  [12],  we  consider 
the  problem 

2 2 2 
- D(x  Du)  + 4x  u ■ - 20x 


Du(0)  - u(l)  - 0, 


which  has  the  solution  u(x) 


5 slnh  2x 
X slnh  2 


The  RRG  approximations  to  u from  several  of  the  Sq(a,v)  subspaces 

were  computed,  and  the  error  graphed  below.  All  computations  were  per- 
formed In  double  precision  on  a POP-10  (with  54  binary  digits).  The 
integrals  required  were  computed  using  composite  Gaussian  quadrature 

with  k-l  nodes  In  each  interval  of  the  mesh.  We  give  the  norms  II  e II  _ 

B 

and  II  De  ||  „ of  the  error  and  its  derivative  (computed  with  composite  k+1 

D 

node  Gaussian  quadrature  rules.) 


CO  LIHEOP  - SLOPE  • -1  99  1 = CO  LINEAP  - SLOPE 
Cl  OUAtPATICS  - SLOPE  ■ -S  28  2 = Cl  OUhC’PATICS  - SLOPE 
C2  cue  ICS  - SLOPE  . -S  89  3 « C2  CUB  ICS  - SLOPE 
C3  OUhPTICS  - SLOPE  • -4  79  4 = C3  OUAPTICS  - SLOPE 
C4  OUIUTICS  - SLOPE  = -5  85  5 » C4  PUINTICS  - SLOPE 


R«TES  OF  CONVERGENCE  --  error  RATES  OF  CONVERGENCE 
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Jesperson  has  solved  Che  quaslllnear  problem 
- D(x  Du)  - - -1^  e“ 

Du(0)  ■ u(l)  “ 0, 

which  has  the  unique  solution  u(x)  - 2 ln( — ^ — ■^)  . Our  theoretical 

8 - X 

results  can  be  made  to  apply  to  this  problem  In  a routine  manner;  cf. 
[14].  The  RRG  approximations  to  u were  computed  using  Newton's  method 
with  an  Initial  guess  of  0;  the  resulting  sequence  of  linear  problems 
was  solved  In  Che  manner  discussed  above.  The  Iteration  was  stopped 

when  the  residual  reached  *^10  In  each  case,  this  required  4 itera- 

tions. 


As  predicted  by  the  theory,  Che  rate  of  convergence  appears  to  be 
k k-1 

h (for  the  error)  and  h (for  the  derivative).  Apparently,  k-1 
quadrature  nodes  per  Interval  are  sufficient  to  obtain  the  predicted 
convergence  rates. 


\ 
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